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Abstract 

A partial differential eigenvalue equation for the density displacement fields 
associated with electronic excitations is derived in the framework of density 
functional theory. Our quantum fluid-dynamical approach is based on a vari- 
ational principle and the Kohn-Sham ground-state energy functional, using 
only the occupied Kohn-Sham orbitals. It allows for an intuitive interpretation 
of electronic excitations in terms of intrinsic local currents that obey a continu- 
ity equation. We demonstrate the capabilities of this non-empirical approach 
by calculating the photoabsorption spectra of small sodium clusters. The 
quantitative agreement between theoretical and experimental spectra shows 
that even for the smallest clusters, the resonances observed experimentally at 
low temperatures can be interpreted in terms of density vibrations. 
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I. INTRODUCTION 



Since its formal foundation as a theory of ground-state properties jl], density functional 
theory has developed into one of the most successful methods of modern many-body theory, 
today also with well-established extensions such as, e.g., time-dependent || and current [JM 
density-functional theory (DFT). In particular in the field of metal cluster physics, DFT 
calculations have contributed substantially to a qualitative and quantitative understanding 
of both ground and excited state properties |||6|]. Understanding the properties of small 
metal particles in turn offers technological opportunities, e.g., to better control catalysis 
0, as well as fundamental insights into how matter grows ||f|. Since the electronic and 
geometric structure of metal particles consisting of only a few atoms still cannot be measured 
directly, photoabsorption spectra are their most accurate probes. Especially the spectra 
of charged sodium clusters have been measured with high accuracy for a broad range of 



cluster sizes and temperatures [K|. A distinct feature of these spectra is that at elevated 



temperatures of several hundred K, in particular for the larger clusters, only a few broad 
peaks are observed, whereas at lower temperatures (100 K and less), a greater number of 
sharp lines can be resolved for clusters with only a few atoms. The peaks observed in 
the high-temperature experiments found an early and intuitive explanation as collective 
excitations in analogy to the bulk plasmon and the giant resonances in nuclei: different 
peaks in the spectrum were understood as belonging to the different spatial directions of the 
collective motion of the valence electrons with respect to the inert ionic background. On the 
other hand, the sharp lines observed in the low-temperature experiments were interpreted 
as a hallmark of the molecule-like properties of the small clusters explicable, in the language 
of quantum chemistry fill , only in terms of transitions between molecular states. 

In this work we present a density functional approach to the calculation of excitations 
that leads us to a unified and transparent physical understanding of the photoabsorption 
spectra of sodium clusters. We first derive a general variational principle for the exact energy 
spectrum of an interacting many-body system. From this, we derive an approximate solution 
in the form of quantum fluid-dynamical differential equations for the density displacement 
fields associated with the electronic vibrations around their ground state. By solving these 
equations, we obtain the eigenmodes within the DFT; hereby only the ground-state energy 
functional and the occupied Kohn-Sham orbitals are required. We demonstrate the accuracy 
of our approach by calculating the photoabsorption spectra of small sodium clusters and 
comparing our results to low-temperature experiments and to configuration-interaction (CI) 
calculations. In this way we can show that also the spectra of the smallest clusters can be 
understood, without knowledge of the molecular many-body wavefunction, in an intuitive 
picture of oscillations of the valence-electron density against the ionic background. 



II. A VARIATIONAL PRINCIPLE 



Starting point for the derivation of the variational principle is the well-known fact that 
for a many-body system described by a Hamiltonian H with ground state |0) and energy E Q , 
the creation and annihilation operators of all the eigenstates obey the so-called equations of 



motion for excitation operators [12 
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(0\O u {H,Oi]\0) = Hw u (0\O u Ol\0) (1) 

(oia[#,a]|o> = ^(oiaaio) = o, (2) 

where O u and 0\, are defined by 

0t| O ) = ]!,), O>> = |0>, and a|0) = 0. (3) 

Of course, the exact solution of these equations are in general unknown. But a variety of 
approximations to the true excited states can be derived from them, e.g., the Tam-Dancoff 
scheme and the small amplitude limit of time-dependent Hartree-Fock theory (RPA). As 



discussed in ||12|| , also higher-order approximations can be obtained. 

Related to these equations, we derive the following variational principle: solving the 
equations ([[]) and (0) for the lowest excited state can be achieved by solving the variational 
equation 

<m. 

in the space of all Hermitean operators. Here, E 3 is defined by 



mi and m 3 are the multiple commutators 

m 1 [Q] = ±(0\[Q,[H,Q}]\0) (6) 
m 3 [Q] = ±(0\[[H,Q],[[H,Q],H]]\Q), (7) 

and Q is some general Hermitean operator that, as will be shown in the course of the 
argument [see Eq. (|15|) below], takes the interpretation of a generalized coordinate. The 
minimum energy E 3 after the variation gives the first excitation energy The second 

excitation with energy hio^ can be obtained from variation in an operator space which has 
been orthogonalized to the minimum Q, and in this way the whole spectrum huj u can be 
calculated. 

The variation 5Q of an operator can be understood as a variation of the matrix elements 
of the operator in the matrix mechanics picture. Therefore, 

= J_ f m 3 [Q] \ * = 1 f m 3 [Q] \ "* J_ f m 3 [Q] \ 

5Q \rnMJ 2 {m^Q]) 5Q \m.lQ]) ' 1 ' 



and noting that the first factors in the expression to the right are just l/(2E; 
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8 f rn 3 [Q] \ 1 6m 3 [Q] m 3 [Q] Sm^Q] 
8Q\m x [Q]) mi[Q] 5Q mi [Q] 2 5Q 



is obtained. With the definition E 3 = hwi, Eq. (§) turns into 
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(10) 



The variations 

5m 3 [Q] = m 3 [Q + SQ] - m 3 [Q] 

5mx[Q] = mi [Q + SQ] - m^Q] (11) 

are evaluated by straightforward application of the commutation rules (||) and (0), leading 
to 

(0\[[6Q,H\, {[H, [H,Q]] - (ftwrfQ) ]|0> = 0. (12) 

With SQ Hermitean, [SQ, H] is anti-Hermitean, and ([12]) therefore is an equation of the form 
c + c* = with 

c=(0\[5Q,H] ([H,[H,Q]]-(hj 1 ) 2 Q)\0)eC. (13) 

Since |0) by definition is the exact ground state of H, and (|13|) must hold for any SQ, the 
equation 

([H,[H,Q)]-(fku 1 ) 2 Q)\0)=0 (14) 

is obtained. It resembles the equation of motion for a harmonic oscillator. Therefore, Q is 
interpreted as a generalized coordinate, and in analogy to the well-known algebraic way of 
solving the harmonic oscillator problem, Q is written as a linear combination 

Q oc 0\ + O x (15) 



of the creation and annihilation operator for the first excited state. Inserting (|T|) into ( 14]) 
leads to the two equations 

[H,[H,O\]]\0) = (hu; 1 ) 2 Ol\0) (16) 

[/f,[ J ff,O 1 ]]|0) = (^ 1 ) 2 O 1 |0) = 0. (17) 

First consider QT5|). After closing with state (1|, one exploits that, by definition, |0) and |1) 
are eigenstates of H and evaluates the outer commutator by letting H act once to the left 
and once to the right. Recalling that (1| = (0|C?i, one finally obtains 

(0\O 1 [HM]\0) = toi(a\OiO\\0). (18) 

This is exactly equation (|l|) for the first excited state. In the same way, @ is obtained from 
(PI), which completes the derivation of the variational principle. 

We would like to point out that in earlier work , the RPA equations have been derived 
with a related technique that made use of both generalized coordinate and momentum 
operators. The advantage of our present derivation is that - although within linear response 
theory - it goes beyond RPA and, due to the formulation in terms of a generalized coordinate 
only, is particularly suitable for the formulation of the variational principle in the framework 
of density functional theory as shown below. 



4 



III. QUANTUM FLUID DYNAMICS FROM THE GROUND-STATE ENERGY 
FUNCTIONAL: A LOCAL CURRENT APPROXIMATION 



In principle, the exact eigenenergies are defined, via Eqs. ([[]), (0), by the variational 
equation ([|), provided that the operator Q is chosen in a sufficiently general form. However, 
just as in the equations of motion technique, one is forced to make some explicit ansatz for 



the form of Q, which necessarily introduces approximations. In Ref. |13| it was shown that 
if Q is taken to be a one-particle-one-hole excitation operator, Eq. (|]) leads to the RPA 
equations. Simplifications of the RPA, in which Q was chosen from restricted sets of local 



operators Q n (r) } were proposed in connection with both semiclassical |H| and Kohn-Sham 



density functionals fl3| . In the present paper, we derive a set of quantum fluid-dynamical 



equations from the variational principle (Q) by choosing Q to a general local operator Q(r). 
These equations are then solved without any restriction other than Eq. ( p3|) below. 

First we recall a relation that is well known in nuclear physics ||15|| : the commutator of 
Eq. ([?[) can be exactly obtained from 

where |a) is the state that results from the unitary transformation 

|a) = e- a5 |0), (20) 

with a being a real and possibly time-dependent parameter, and S the so called scaling 
operator defined by 

S=[H,Q}. (21) 

Assuming that Q is just a function of r and that the potentials in H do not contain derivatives 
with respect to r, as is the case for Coulombic systems, Eq. ( |2~T| ) is easily evaluated: 

N a A^c 

2 



N c N e 

S = J2 = E 9 ( V * U W) + U W ' V *- ( 22 ) 



i=l i=l 

Here, the displacement field 



u(r) = --VQ(r) (23) 



m 



has been introduced, and N e is the number of electrons. 

These equations can be related to DFT by noting that, first, we can introduce a set of 
single particle orbitals {^(i - ;)}, and from the scaled single particle orbitals, a scaled single 
particle density can be constructed via 



N L 

|2 



n{r, a ) 



J2 |e~ Qs(r) Vv(r)| = e- aSn n(r), (24) 



with a density scaling operator 
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S n = (Vu(r)J + u(r) • V. (25) 
Second, Eq. (^) can straightforwardly be evaluated for a local Q(r), 

mi[Q] = ^ / u(r) • u(r)n(r) d 3 r, (26) 

showing that m x depends only on n and u and is similar to a fluid-dynamical inertial 
parameter. And third, we replace the expectation value in Eq. (|19|) by 

m 3 lQ] = ~(a\H\a)l^ 1 -^EHr,a)] (27) 



a=0 



where E[n] is the usual ground-state Kohn-Sham energy functional 

E[n; {R}] = T s [n] + £ xc [™] +jj J "y^T <*V ^ + / n(r)V ion (r; {R}) d 3 r. (28) 

Eq. (|26| ) is exact and also Eq. (|2H]) can be verified order by order, but Eq. (p7[) goes beyond 
the safe grounds on which the energy functional is defined. However, the replacement of 
an energy expectation value by the energy functional is intuitively very plausible, and its 
practical validity can be judged a posteriori by the results. A further strong argument for 
why really the density should be the basic variable can be made by calculating the derivative 
with respect to time of the scaled density, using Eqs. ([24]) and fl25|). 



^-n{r, a(t)) = -S n a(t) n{r, a(t)) = -V[d(t) u(r) n{r, a(t))), (29) 

where for the sake of clarity we now explicitly wrote the time dependence of a. Since 

j(r,*) = d(*)u(r)n(r,a(*)), (30) 

is a current density, Eq. ((29^) is just the continuity equation dn(r, a(t))/dt + Vj(r,t) = 0. 
Thus, the variational principle Eq. @ with a local function Q(r) describes excitations by 
intrinsic local currents. The time dependence of the parameter a is obviously harmonic, i.e., 
a(t) oc cos(tovt), since the present derivation is based on linear response theory. 

The physical significance of the variational approach now being clear, it remains to derive 
the actual equations that determine the displacement fields u(r) and the energies hw that 
are associated with particular excitations. Starting from Eq. (|10"D and using an explicit 
notation, 

5m 3 [u[Q{r)]\ ^^ 2 <Wu[Q(r)]] = Q= r^, pm 3 [n(r)] _ ^^ 2 5m x [u(r)] \ 6u{r") 



5Q(r') v 17 5Q(r>) J { 5u(r") v LJ 5u{r") J 5Q(r') 

(31) 

follows by virtue of the chain rule for functional derivatives. Thus, solutions of 

zgm = (^ s -^i (32) 

Our Our 
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will also be solutions to Eq. (|I~0"D and thus Eq. (|j). ni\ is already given as the functional 
] by Eq. (H), and 777,3 [u] is readily obtained by inserting the scaled Kohn-Sham orbitals 
and density from Eq. ( p4|) into the energy functional Eq. ( p8|) and calculating the second 
derivative with respect to the parameter a, Eq. (|27|) . The final equations are then derived in 
a lengthy but straightforward calculation from Eq. (^) by explicitly performing the variation 
on u. Using the usual definition 



5ms[u(r)] <5m 3 [u](r) 5m 3 [u(r)] 5m 3 [u(r)] 
5u(r') ~ 



8uJr' 



-e x + 



5u y (r' 



5uJr' 



(33) 



where are the unit vectors in the Cartesian directions, a set of three coupled, partial differ- 
ential eigenvalue equations of fourth order for the Cartesian components Uj(r) is obtained: 



5m 3 [u] 
5uj(r) 



2 5mi[u] 
8uj{r) ' 



J = 1,2,3, 



(34) 



where 



<5mi[u] m 

— ^ = ¥ n(r) Uj (r), 



(35) 



<5m 3 [u] £777,3 [u] 5mf [u] 5rrq 2 [u] Sm^ [u] 
5uj{v) 5uj{r) 8uj{r) 8uj{r) Suj(r) 



(36) 



and 



5m| in [u] 



~2m2 



N c 3 



( A ^mj (d j u i )(d i il)^ n ) + (djdiUi)il)* m + u^djdiipl 

m=l i=l *■ 

(djUi)(diAi/; m ) + u^djdiA^m) tp* m - ui 



+ 



+2 



(djipm) A[-(diUi)ip m + Ui(diip ri 



djA[-(diUi)ip m + Ui(diip r , 



(37) 



5mf [u 
Suj{r) 



1 3 

IE 



i=l 



n((9 i 'U l )(9ifKs) - {diUi){djV K s) J +Wt(^(^9,f KS ) - (din)(djV KS 



(38) 



5w,-(r) 



n 



^(^(r'))n(r') + n,(r')(^(r')) 



rj rj dV, 



i=l 



r — r 



/|3 



(39) 



5mg c2 [u] 
Suj(r) 



-n^2 (dj({diUi)n + Ui(din))^^- + ({diu^n + w^n)) (^^^) 



(40) 
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where we used the shorthand notation d\ = d/dx etc., and indicated the terms to which 
derivatives refer by including them in parenthesis. The usual Kohn-Sham and exchange- 
correlation potential are denoted by t>Ks an d ^xc, respectively. 

Eqs. ( 53]) - ( f40|) are our quantum fluid-dynamical equations. In analogy to the local 
density approximation (LDA) used for v xc , we term our scheme the local current approxi- 
mation (LCA) to the dynamics, due to the use of a local function Q(r) in the variational 
principle (|). It should be noted that the above equations differ from the equations derived 
earlier in a semiclassical approximation [13] or by explicit particle-hole averaging |I3|. Due 



to the fact that our approach is completely based on the Kohn-Sham density functional and 
therefore contains the full quantum-mechanical shell effects in the ground-state density, it is 
also different from some fluid-dynamical approaches developed in nuclear physics [|16j (and 
used in cluster physics fl?]]) which involved either schematic liquid-drop model densities or 
semiclassical densities derived from an extended Thomas-Fermi model. 

Although Eqs. ([34j) - ([!(]) look rather formidable, they can be solved numerically with 
reasonable computational effort, and we have done so for the sodium clusters Na2 and Najj~. 
The Kohn-Sham equations were solved basis-set free on a three-dimensional Cartesian real- 
space grid using the damped gradient iteration with multigrid relaxation The ionic 
coordinates were obtained by minimizing the total energy using a smooth-core pseudopo- 
tential 0. For E xc , we employed the LDA functional of Ref. [[l!|. The Uj(r) were expanded 
in harmonic oscillator wavefunctions and we explicitly enforced Eq. (^). The convergence 
rate of the expansion can be improved by adding a few polynomial functions to the basis. 
By multiplying Eqs. ( [32]) and subsequently (]34])-(f40|) from the left with u and integrating 
over all space, a matrix equation for the expansion coefficients is obtained which can be 
solved using library routines. The square roots of the eigenvalues then give the excitation 
energies and from the eigenvectors, the oscillator strengths can be computed. 

It should be noted that for systems as small as the ones studied here, generalized gra- 
dient approximations |2(J and their extensions ^TJ in general are a better approximation 
to the exchange and correlation energy than the LDA. However, in the present case LDA 
is a good approximation since the valence electrons in sodium clusters are strongly delocal- 
ized. Furthermore, the bond length underestimation in LDA which was shown to strongly 
influence optical properties p2|p5[ is corrected to a large extent by using the smooth-core 
pseudopotential which by construction gives bondlegths close to the experimetal ones when 
used with the LDA §. 

Fig. [I] shows the experimental photoabsorption spectrum [|23J of Na2 in the upper left 
picture (adapted from Ref. 0), and below the spectrum obtained in the just described LCA. 
We introduced a phenomenological line broadening in the LCA results to guide the eye. The 
LCA correctly reproduces the electronic transitions, despite the fact that only two electrons 
are involved. Due to Eq. fl2~9D , one can very easily visualize how the electrons move in a 
particular excitation by plotting the corresponding Vj(r), giving a "snapshot" picture of 
dn/dt. For the two main excitations, a crossection of this quantity along the symmetry axis 
(z axis) is shown in the lower left and upper right contourplots, and the ground-state valence 
electron density is shown in the lower right for reference. In the plots of dn/dt, shadings 
darker than the background grey indicate a density increase, lighter shadings indicate a 
decrease. It becomes clear that the lower excitation corresponds to a density oscillation 
along the z axis whereas the higher excitation corresponds to two energetically degenerate 



S 



oscillations perpendicular to the symmetry axis. (For the sake of clarity, we plotted the 
corresponding oscillator strengths on top of each other in the photoabsorption spectrum.) 
This is exactly what one would have expected intuitively. But the plots reveal that besides 
the expected general charge transfer from one end of the cluster to the other, the presence 
of the ionic cores hinders the valence electrons to be shifted freely, creating a density shift 
of reverse sign in between the ionic cores. 

Fig. |2] shows the ionic ground-state configuration of Na^ with our labeling of axes in the 
upper left, the experimental low-temperature (~ 100 K) photoabsorption spectrum |10| in 



the upper right, the LCA photoabsorption spectrum in the lower left, and the CI spectrum 



adapted from Ref. [|TTJ in the lower right. Again, a phenomenological line broadening was 
introduced in the presentation of both the LCA and the CI results. The LCA spectrum 
again is in close agreement with the experimentally observed spectrum, showing three intense 
transitions. With our choice of the coordinate system, the lowest excitation corresponds to 
a density oscillation in z direction, whereas the two higher excitations oscillate in both x 
and y directions. In the interpretation of the LCA results, it must be kept in mind that 
due to our finite grid spacing the numerical accuracy for the excitation energies is about 
0.03 eV, which is absolutely sufficient in the light of the physical approximations that we 
are making. But due to this finite numerical resolution and the fact that we evaluate each 
direction of oscillation separately, the x and y components of the excitations at 2.7 eV and 
3.4 eV, which really should be degenerate for symmetry reasons, appear as extremely close- 
lying double lines. However, since the symmetry of the cluster was in no way an input to our 
calculation, it is a reassuring test that the LCA, indeed, fulfills the symmetry requirement 
within the numerical accuracy. Furthermore, it is reassuring to see that with respect to 
the relative heights of the peaks the LCA is very close to the CI results, with differences 
observed only in the small subpeaks that are not seen experimentally anyway. And small 
differences to the CI calculation are already to be expected simply because of the use of 
different pseudopotentials and the resulting small differences in the ground-state structure. 



IV. CONCLUSION 



In summary, we have derived a set of quantum fluid-dynamical equations from a general 
variational principle for the excitations of a many-body system. The equations describe 
here the eigenmodes of the system's (valence) electrons and require only the knowledge of 
the occupied ground-state Kohn-Sham orbitals. From these equations, we have computed 
the photoabsorption spectra for small sodium clusters and find quantitative agreement with 
the experimentally observed peak positions. Thus, even low-temperature photoabsorption 
spectra can be understood in an intuitive picture of density oscillations, without knowledge 
of the true (or any approximate) many-body wavefunction. 
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FIGURES 




FIG. 1. From top left to bottom right: Experimental photoabsorption cross section (241] and 
LCA cross section S of Na2 in arbitrary units versus eV, contour plots of the density change 
associated with the first excitation, the second excitation, and contour plot of the ground-state 
valence electron density. The length unit for the axes of the contour plots is the Bohr radius qq. 
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FIG. 2. Upper left: ionic ground-state configuration of Na^", lower left: corresponding LCA 
photoabsorption spectrum, upper right: experimental low-temperature photoabsorption spectrum 
PI , lower right: Configuration-Interaction photoabsorption spectrum from Ref. [11|. See text for 
discussion. 
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